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1  Introduction 


This  report  documents  the  research  conducted  under  AFOSR  Grant  F49620-92- 
J-0488,  “Compressible  Flow  Turbulence  Simulation  and  Modeling  via  Additive 
Turbulent  Decomposition,”  during  the  period  09/01/92-11/30/95. 

We  begin  by  noting  that  four  graduate  students  have  been  supported  during 
the  course  of  this  grant:  Mr.  D.  Denger,  Ms.  S.  Flynn,  Mr.  E.  C.  Hylin  and 
Mr.  D.  C.  Weatherly.  The  first  two  chose  to  leave  the  project  before  the  end 
of  the  first  year;  Ms.  Flynn  received  her  M.S.  in  Mechanical  Engineering  shortly 
thereafter,  and  Mr.  Denger  is  expected  to  receive  his  Ph.D.  in  Mechanical  Engi¬ 
neering  at  least  by  the  end  of  Summer  1996.  Both  of  the  remaining  two  students, 
Mr.  Hylin  and  Mr.  Weatherly,  anticipate  completing  their  Ph.D.s  in  Mechanical 
Engineering  during  the  current  semester.  Spring  ’96.  Thus,  the  grant  will  have 
aided,  very  crucially  in  the  latter  two  cases,  in  producing  four  advanced  engineer¬ 
ing  degrees. 

The  work  to  be  reported  here  deviates  somewhat  from  that  originally  proposed, 
at  least  in  part  due  to  our  evolving  understanding  of  the  additive  turbulent  de¬ 
composition  (ATD)  formalism.  In  particular,  it  was  originally  proposed  to  study 
both  the  simulation  and  modeling  aspects  of  ATD  applied  to  compressible  flow. 
However,  by  the  beginning  of  the  second  year  of  the  project,  parallelization  stud¬ 
ies  of  ATD  being  conducted  under  NSF  sponsorship  were  beginning  to  show  that 
while  predicted  parallel  speedups  were  indeed  achievable,  the  number  of  proces¬ 
sors  needed  for  high-Re  engineering  calculations  far  exceeded  that  likely  to  be 
available  in  the  near  term,  so  the  direct  simulation  version  of  ATD  is  not  likely 
to  be  widely  useful  very  soon.  At  the  same  time,  and  in  part  motivated  by  these 
findings,  we  began  incompressible  flow  studies  employing  ATD  as  the  foundation 
for  a  new  class  of  turbulence  modeling  techniques.  Initial  results  were  extremely 
promising,  so  we  pursued  this  approach  for  both  compressible  and  incompressible 
flow. 

The  basic  idea  underlying  this  approach  has  been  reported  in  earlier  progress 
reports  for  this  and  other  grants,  and  it  will  be  discussed  at  some  length  in  the 
sequel.  It  consists  of  solving  unaveraged  equations  for  the  large-scale  quantities — 
these  equations  containing  both  large-  and  small-  (and  any  appropriate  inter¬ 
mediate-)  scale  variables — and  modeling  everything  else  using  nonlinear  algebraic 
maps.  This  approach  is  extremely  straightforward.  It  greatly  reduces  the  num¬ 
ber  of  quantities  that  must  be  modeled,  compared  with  Reynolds-averaged  ap¬ 
proaches,  and  it  allows  time-accurate  calculations  including  interactions  between 
turbulence  and  other  physical  phenomena. 

It  is  crucial  to  recognize  that  it  is  the  fluctuating  primitive  variables  themselves, 
and  not  their  correlations,  that  are  being  modeled.  Figure  1  depicts  some  typical 
results  for  2-D  incompressible  flow  over  a  backward-facing  step  at  Re  =  1.32  x  10®. 
Parts  (a)  and  (b)  show  instantaneous  streamlines  of  the  complete  and  small-scale 
velocity  fields  respectively.  Of  particular  interest  in  part  (a)  of  the  figure  are  the 
secondary  vortices  at  the  downstream  end  of  the  main  vortex.  It  is  important  to 
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Figure  1:  Modeled  flow  over  a  2-D  backward-facing  step. 


note  that  these  are  transient,  and  that  they  do  not  occur  if  the  small-scale  model 
is  omitted.  Part  (b)  of  the  figure  indicates  that  the  instantaneous  maiss-conserved 
small-scale  velocity  exhibits  a  wide  range  of  length  scales,  at  least  some  of  which 
interact  with  the  large-scale  motion. 

Finally,  in  part  (c)  we  display  the  time-averaged  behavior.  The  salient  points  to 
note  are:  i)  the  main  vortex  is  smooth  and  of  a  length  comparable  to  experimental 
observations,  and  ii)  a  secondary  vortex  in  the  lower  corner  of  the  step  has  been 
resolved,  even  on  the  fairly  coarse  (but  stretched)  61  x  46  grid.  It  is  also  worth 
mentioning  that  Reynolds  stresses  and  turbulent  kinetic  energy  are  in  reasonable 
agreement  with  experimental  data. 

The  run  time  required  for  this  simulation  (without  parallelization)  is  within 
a  factor  of  1.5  of  that  required  by  the  standard  k-t  model  on  the  same  grid. 
Fairly  straightforward  parallelization  would  reduce  this  considerably,  even  if  only 
relatively  few  processors  are  employed  (as  few  as  two). 

In  the  following  sections  of  this  report  we  will  present  the  chaotic-map  mod¬ 
eling  procedure  in  fair  detail  (a  complete  treatment  is  available  in  Hylin  &  Mc¬ 
Donough  [2]),  along  with  discussions  concerning  2-D  and  3-D  implementations 
for  the  GASP  compressible  code.  Computed  results  for  a  Mo©  =  2.85  flow  over 
a  compression  ramp  will  be  presented  for  the  2-D  case,  and  initial  3-D  results 
associated  with  vortex-stimulated  flows  will  also  be  given.  Completion  of  the  3-D 
incompressible  and  3-D  compressible  versions  of  the  ATD/chaotic-map  modeling 
formalisms  is  the  remaining  Ph.D.  dissertation  research  requirement  for  Mr.  Hylin 
and  Mr.  Weatherly  respectively. 

2  Review 

Since  September  of  1994,  we  have  continued  to  refine  our  idea  of  a  new  class  of 
turbulence  models,  which  model  fluctuations  in  the  primitive  variables  instead 
of  modeling  statistically  averaged  quantities  such  as  the  Reynolds  stress  or  the 
turbulent  kinetic  energy  and  dissipation.  By  taking  this  approach,  the  subgrid- 
scale  stresses  generated  by  our  models  are  guaranteed  to  be  physically  realizable — 
unlike  many  classical  turbulence  models,  ours  will  never  generate  a  negative  kinetic 
energy.  In  fact,  it  has  proven  possible  to  construct  models  that  by  design  duplicate 
many  features  of  physical  turbulence  at  small  scales. 

Our  models  are  enabled  by  several  fundamental  ideas.  First,  to  decompose  the 
conserved  variables  in  the  equations  of  fluid  motion  into  large-scale  and  small-scale 
components,  and  then  to  decompose  the  equations  themselves  into  two  coupled 
sets  of  equations  governing  the  evolution  of  these  components.  Second,  to  replace 
the  evolution  equations  for  the  small-scale  quantities  with  stochastic  small-scale 
fields;  that  is,  to  create  stochastic  models  for  the  small-scale  fluctuations  in  the 
conserved  variables,  and  use  these  models  to  generate  the  spatio-temporal  fields  as¬ 
sociated  with  these  quantities — for  example,  the  small-scale  velocity  field.  Third, 
to  generate  the  necessary  stochastic  variables  via  chaotic  maps,  and  to  tailor  these 
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so  as  to  mimic  as  closely  as  possible  the  dynamical  and  statistical  behavior  of  a 
physical  small-scale  flow.  Fourth,  to  scale  the  amplitude  of  the  small-scale  fluctu¬ 
ations  to  match  the  turbulent  spectra,  especially  the  spectrum  for  the  turbulent 
kinetic  energy.  And  fifth,  to  use  local  samples  of  the  large-scale  flow  structure  to 
estimate  the  local  anisotropy  and  one-point  velocity  correlations  for  the  subgrid- 
scale  flow. 

In  our  previous  report,  certain  areas  of  our  models  were  identified  as  being 
most  in  need  of  refinement.  These  included  the  decomposition  of  the  governing 
equations,  the  relation  between  large-scale  flow  parameters  and,  the  dynamical 
parameter(s)  governing  our  chaotic  maps,  the  scaling  in  compressible  flows  of  the 
energy  and  density  fluctuations,  and  the  form  of  the  chaotic  maps  as  it  relates  to 
the  statistical  properties  of  real  turbulence.  We  have  made  significant  progress 
in  two  of  these  areas:  the  decomposition  of  the  equations  and  the  form  of  the 
chaotic  maps.  Refinement  of  the  energy  and  density  fluctuations  still  awaits  better 
turbulent  spectra  for  these  quantities,  while  an  improved  relation  between  the  flow 
and  map  parameters  likely  depends  on  progress  in  dynamical  systems  theory. 

3  Decomposition 

Previously  we  had  split  the  large-scale  and  small-scale  equations  on  something  of 
an  ad  hoc  basis,  relying  on  heuristic  arguments  that  the  large-scale  and  small- 
scale  equations  should  have  the  same  form,  and  that  they  should  take  the  form  of 
transport  equations  for  the  large-scale  and  small-scale  quantities.  There  was  some 
evidence,  however,  that  the  splitting  could  be  made  more  mathematically  rigorous 
if  it  was  defined  in  terms  of  projections,  such  that  the  large-scale  equation  was 
derived  via  a  projection  onto  the  large-scale  basis,  and  the  small-scale  equation 
derived  via  a  complementary  projection.  We  have  succeeded  in  the  effort  to  define 
such  a  splitting,  which  will  now  be  described. 

3.1  Projection  Operators  for  the  Large  and  Small  Scales 

In  additive  turbulent  decomposition,  as  in  large-eddy  simulation,  the  first  step  is  to 
split  the  dependent  variables  into  large-scale  and  small-scale  components  by  means 
of  a  spatial  filter.  If  for  some  set  of  basis  functions  this  filter  is  spectrally  sharp, 
then  the  large-scale  components  represent  the  first  few  terms  in  the  corresponding 
functional  expansion  for  the  variables,  while  the  small-scale  components  represent 
the  remainder  of  the  expansion.  Thus  if  the  velocity  u  may  be  expressed  as  the 
weighted  sum  of  a  series  of  basis  functions: 

«  =  (1) 

K 

where  4>{k‘x)  is  the  basis  function  for  the  wavenumber  k,  operating  on  the  position 
vector  X,  while  uk,  is  the  corresponding  weight;  then  the  large-scale  velocity  may 


4 


be  expressed  as 


=  5^  UK<f>{K-x)  -  Pl(w),  (2) 

|«|<Ar 

while  the  small-scale  velocity  is  expressed  as 

Us=  UK<f>{K-x)  =  {I  -rL){u)  =Fs{u).  (3) 

\K\>N 

Equations  (2),and  (3)  define  Pf,  as  the  projection  operator  for  the  large-scale  basis, 
and  P5  as  the  complementary  projection  operator  for  the  small-scale  basis. 


3.2  The  Large-  and  Small-Scale  Equations 

If  the  projection  operators  Pl  and  P5  are  applied  to  the  incompressible  Navier- 
Stokes  equations,  then  since  the  time  derivative  operator  commutes  with  the  spa¬ 
tial  projectors,  one  obtains  differential  equations  for  the  large-scale  and  small-scale 
velocities: 

+  Pl  -I-  V‘{usUl)  +  V'{ulus)  -f  V-(w5Us)] 

=  — Pz,(VP) -1- t'P£,A(ui, -t- 1*5),  (4) 

dus 

-f  P5[V*(rtx,«i)  -1-  V-{usUl)  +  V‘(ulUs)  -f-  V*(W5U5)] 

=  — P5(VP) -f  1/P5A(w2, -f  1*5),  (5) 

Pl[V»(wl -I- 145)]  =  0,  (6) 

and 

P5  [V*(w£,  +  W5)]  =  0.  (7) 

Defined  in  this  way,  this  coupled  system  of  equations  for  ul  and  Us  is  equivalent 
to  the  Navier-Stokes  equations.  Because  the  projection  and  convolution  operators 
do  not  commute,  there  is  necessarily  coupling  through  the  advective  terms,  but 
in  general  there  is  coupling  through  the  diffusive  terms  as  well. 

An  exception  occurs  when  the  Fourier  functions  are  used  for  the  basis  func¬ 
tions  (f).  Then  the  projections  commute  with  spatial  differentiation;  the  diffusive 
terms  uncouple,  as  do  (6)  and  (7),  so  the  incompressible  continuity  equation  holds 
separately  on  the  large  and  small  scales.  Also  in  this  case  a  large-scale  and  a  small- 
scale  pressure  may  be  defined  straightforwardly  through  corresponding  pressure 
Poisson  equations: 

~  +  ^LiusUi)  +  PL(wz,tt5)  +  PL(w5tt5)]}  ,  (8) 
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and 


^Ps  —  — v«  {V*  [Ps(WLUi)  +  Fs{usUl)  +  Ps(mi,'U5)  +  P5(Msti5)]}  •  (9) 

Certain  other  bases — the  Chebyshev  and  Legendre  polynomials,  for  example — 
permit  a  partial  simplification  of  the  large-scale  and  small-scale  equations.  How¬ 
ever  the  small-scale  velocity  always  influences  the  large-scale  advection,  as  the 
large-scale  velocity  influences  the  small-scale  advection,  and  consequently  the 
large-scale  and  small-scale  velocities  must  (in  principle)  be  solved  for  simulta¬ 
neously. 

4  Stochastic  Small-Scale  Fields 

Before  proceeding  with  the  description  of  our  subgrid-scale  model,  let  us  consider 
some  features  of  a  stochastic  small-scale  field.  If  the  large-scale  velocity  field  is 
defined,  for  example,  on  a  finite-difference  grid  with  a  spacing  of  h  (the  coarse 
grid),  the  small-scale  field  may  be  defined  on  a  grid  with  a  spacing  of  /i/2  (the 
fine  grid).  Since  in  this  case  the  fine  grid  includes  the  points  of  the  coarse  grid, 
the  basis  defined  by  the  fine  grid  overlaps — indeed,  it  includes — that  defined  by 
the  coarse  grid.  However,  what  distinguishes  the  velocity  fields  we  shall  define  on 
the  fine  grid  as  being  stochastic  small-scale  fields  are  the  following  features: 

1.  A  stochastic  variable  or  variables  will  be  used  to  assign  a  value  to  each 
gridpoint. 

2.  These  values  will  be  updated  at  time  intervals  determined  in  part  by  the  grid 
spacing,  and  commensurate  with  the  expected  turnover  time  for  subgrid- 
scale  eddies. 

3.  The  maximum  magnitude  of  these  values  will  be  a  small  fraction  of  the 
maximum  large-scale  velocity. 

4.  Over  distances  larger  than  a  small  multiple  of  the  grid  spacing,  these  values 
will  be  uncorrelated. 

5.  The  large-scale  velocity  field  will  be  filtered  to  reduce  or  eliminate  fluctua¬ 
tions  occuring  on  a  scale  smaller  than  a  small  multiple  of  the  grid  spacing. 

As  a  consequence  of  these  features,  the  large-scale  and  small-scale  velocity  fields 
evolve  on  different  spatial  and  temporal  scales.  And  while  transitory  large-scale 
structures  may  appear  in  the  small-scale  field,  they  will  endure  only  a  small  frac¬ 
tion  of  the  characteristic  large-scale  time,  and  will  have  an  amplitude  that  is  only 
a  small  fraction  of  the  characteristic  large-scale  velocity.  Nonetheless,  it  is  such 
transitory  structures  that  will  provide  the  means  for  the  advective  coupling  of  the 
large  and  small  scales. 
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5  Subgrid-Scale  Model 

The  idea  we  are  following  is  to  create  a  stochastic  small-scale  field  by  defining 
a  fluctuating  value  at  each  point  in,  say,  a  flnite-difference  grid.  (We  have  also 
applied  this  technique  to  a  finite- volume  discretization.)  From  a  mathematical 
viewpoint,  each  of  these  fluctuating  values  can  generally  be  resolved  into  the 
sum  of  several  contributing  terms,  each  of  which  takes  the  form  of  an  amplitude 
coefficient  multiplying  a  normalized  stochastic  variable.  Our  analysis,  however, 
proceeds  from  a  physical  viewpoint  in  which  each  fluctuating  value  is  produced 
by  a  combina4;ion  of  physical  effects,  which  are  added  or  multiplied  together  as 
appropriate.  We  shall  now  consider  these  components  one  at  a  time. 


5.1  Basic  Forms  and  Amplitudes 

For  compressible  flow,  the  model  for  the  subgrid-scale  velocity  fluctuations  Ws 
takes  the  form 


Ws  =  C^w^ReKCo  M).  (10) 

Here  Ws  is  the  subgrid-scale  velocity  weighted  by  the  square-root  of  the  density, 
while  on  the  right-hand  side  the  first  three  terms  together  specify  the  amplitude 
of  the  fluctuations,  the  fourth  term,  is  a  vector  of  anisotropy  coefficients,  and 
the  fifth  term,  AT  is  a  vector  of  values  from  the  chaotic  maps  that  generate  the 
normalized  stochastic  variables.  These  last  two  terms  are  multiplied  together  in 
a  vector  Hadamard  product,  defined  for  two  vectors  and  a  unit  vector  i  according 
to: 


(CoAT).i=(C-i)(AT.i).  (11) 

For  incompressible  flow  velocities  a  similar  expression  is  used,  in  which  a  quantity 
u,,  a  friction  velocity,  replaces  the  quantity  w,,  the  friction  velocity  weighted  by 
the  square-root  of  the  density,  and  the  result  is  a  pure  velocity,  us- 

The  terms  in  the  expression  for  the  velocity  fluctuations  have  the  following 
meanings.  Rch  is  a  small-scale  Reynolds  number. 


Reh  = 


u  ’ 


(12) 


where  h  is  the  local  grid-cell  characteristic  length  (which  in  generalized  coordinates 
we  take  to  be  the  matrix  2-norm  of  the  inverse  of  the  coordinate  transformation 
matrix),  i/  is  the  (molecular)  kinematic  viscosity,  and  1|  Vtt||  denotes  the  2-norm 
of  the  Jacobian  matrix  (with  respect  to  spatial  coordinates)  of  the  large-scale 
velocity  vector.  As  mentioned  above,  iv,  is  a  density-weighted  friction  velocity 
and  is  defined  according  to: 


=  /32(i/||Vw|])2  =  p-iu. 


(13) 
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Cu  is  an  amplitude  pre-factor  that  in  general  is  a  function  of  the  small-scale 
Reynolds  number  Rch  and  the  Taylor  microscale  Reynolds  number  Rex,  but  in 
the  inertial  subrange  assumes  a  nearly  constant  value:  ^  0.877. 

I 

The  amplitude  factor  (CuW^Rel)  determined  for  the  velocity  fluctuations  is  in 
accordance  with  Komogorov’s  second  similarity  hypothesis  [6];  this  hypothesis  im¬ 
plies  that  over  a  distance  h  in  the  inertial  subrange,  turbulent  velocity  fluctuations 
Au  should  scale  as 

Au~(he)3,  (14) 

where  c  is  the  energy  dissipation. 

In  the  compressible  case,  we  compute  additionally  the  density  and  energy 
fluctuations.  The  former  is  given  by 

PS  =  pL CpCu  Rel^Mp.  (1 5) 


Here  Mp  is  an  instance  of  the  normalized  stochastic  variable,  while  the  preceding 
terms  together  constitute  the  amplitude  factor.  Several  of  these  terms  have  been 
previously  defined;  of  the  others,  pi  is  the  large-scale  density,  L  is  a  large-scale 
characteristic  length,  and  Cp  is  a  possible  0(1)  coefficient.  This  amplitude  factor 
was  derived  by  assuming  that  the  ratio  of  the  small-scale  density  fluctuations  to 
the  large-scale  density  is  of  the  same  order  as  the  ratio  of  the  small-scale  velocity 
fluctuations  to  the  large-scale  velocity.  That  is. 


£5  ^  ||«5|| 

PL  IIwlII' 


(16) 


Rather  than  use  such  an  ad  hoc  assumption,  we  would  prefer  to  obtain  the  density 
amplitude  function  in  the  same  way  we  obtained  the  velocity  amplitude  function: 
by  matching  an  assumed  power-law  scaling  function  to  an  experimental  and/or 
theoretical  spectrum.  Unfortunately,  as  far  as  we  know,  such  a  spectrum  has  yet 
to  be  determined.  At  least,  the  assumption  we  have  adopted  does  not  appear  to 
lead  to  any  obvious  contradictions. 

For  the  energy  amplitude  functions,  the  basic  assumption  is  that  the  small- 
scale  internal  and  kinetic  energy  fluctuations  are  of  the  same  order;  that  is. 


2{p^)s 


u  ~  {pe)s  ~  {pE) 


5’ 


(17) 


where  e  is  the  internal  energy  and  E  is  the  total  energy.  This  assumption  is 
supported  by  some  work  by  Kida  &  Orszag  [5].  However,  we  expect  there  to  be  two 
sources  of  small-scale  kinetic  energy  fluctuations:  a  purely  small-scale  term  and  a 
cross-term  involving  both  the  large-scale  and  small-scale  velocities.  Accordingly, 
the  expression  we  assume  to  model  the  small-scale  energy  fluctuations  in  total 
energy  is: 

{pE)s  =  AecMcMs  +  Aes^si  (18) 
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where  Aec  and  Aes  are  the  amplitude  factors  for  the  cross  term  and  purely  small- 
scale  term  respectively,  while  Me  and  Ms  are  two  different  independent  instances 
of  the  normalized  stochastic  variable.  The  amplitude  factors  are  given  by: 

Aes  =  CesCIwIRcI  ,  (19) 

Aec  =  CEcCuwlRel.  (20) 

On  the  basis  of  the  assumption  that  the  fluctuations  in  the  kinetic  and  internal 
energy  should  be  of  the  same  order,  we  expect  that  |  <  Ces  ^  3,  and  Ces  ^ 
Cec  ^  ‘^Ces-  As  in  the  case  of  the  density  fluctuations  we  should  prefer  to  base 
these  amplitudes  on  an  experimental  or  theoretical  spectrum,  but  such  seem  to 
be  lacking.  The  experimental  determination  of  spectra  for  turbulent  density  and 
internal-energy  fluctuations  in  compressible  flow  is  much  needed. 

5.2  Anisotropy 

While  the  amplitude  of  the  small-scale  density  fluctuations  can  (presumably)  be 
expressed  by  a  single  coefficient,  and  the  amplitude  of  the  energy  fluctuations 
by  one  coefficient  on  each  of  the  two  contributing  terms,  the  amplitude  of  the 
small-scale  velocity  fluctuations  must,  in  general,  be  expressed  by  a  vector  of 
coefficients,  because  the  velocity  is  a  vector  function.  We  have  chosen  to  express 
this  amplitude  vector  as  the  product  of  an  amplitude  coefficient  and  a  vector  s  of 
anisotropy  coefficients. 

There  are  two  potential  sources  of  anisotropy,  one  arising  from  the  physics  and 
the  other  from  the  numerics  of  the  flow.  The  former  is  simply  the  anisotropy  of  the 
physical  flow  structure,  which  is  manifest  in  a  difference  between  the  magnitudes 
of  the  three  components  of  the  turbulent  fluctuations  in  the  physical  velocity;  the 
latter  is  due  to  whatever  anisotropy  there  is  in  the  grid  spacing.  Both  impact  the 
amplitudes  that  should  be  specified  for  the  velocity  fluctuations. 

In  modelling  the  anisotropies  inherent  in  the  flow,  we  hypothesize  that,  locally 
in  space,  the  anisotropy  is  a  smooth  function  of  the  scale.  Thus,  at  each  location 
in  space,  the  anisotropy  associated  with  eddies  just  larger  than  the  grid  scale  is 
almost  the  same  cis  the  anisotropy  present  in  the  largest  eddies  on  the  subgrid 
scale.  Considering  the  fluctuations  produced  by  our  model  to  be  equivalent  to  the 
latter,  we  conclude  that  we  may  obtain  the  anisotropy  coefficients  by  sampling 
the  local  structure  of  the  flow  just  above  the  grid  scale.  The  flow  at  this  scale  is 
determined  by  running  the  large-scale  velocity  field  through  a  high-pass  filter. 

It  seems  reasonable  to  base  the  anisotropy  coefficients  on  the  root-mean-square 
values  of  the  local  velocity  fluctuations.  Let  the  filtered  velocity  at  a  point 
(xi,a:2,a:3)  be  given  by  w  =  (ui,  it2,  U3),  and  consider  the  differences  from  this 
velocity  at  points  a  distance  S  =  0{h)  away.  In  particular,  consider  the  velocity 
fluctuation 


Aui  =  ui{x  +  S)  —  ui{x). 


(21) 


9 


It  can  be  shown  that 


{(Au)=)  =  |||Vu,UV  +  0(i*).  (22) 

It  seems  reasonable  to  base  the  anisotropy  coefficients  on  the  root-mean-square 
values  of  the  local  velocity  fluctuations.  Therefore,  using  the  above  estimate,  we 
set 


Sk 


=  V3 


2\2 


((Auf)' 


^/3■ 


liv^ii  • 


(23) 


If  the  flow  is  being  simulated  on  a  cubical  grid,  the  anisotropy  coefficients  are 
completely  specified  by  (23),  but  such  will  rarely  be  the  case.  Nevertheless,  in 
the  case  of  a  potentially  anisotropic  grid,  where  the  flow  is  described  in  terms  of 
generalized  coordinates,  (23)  is  still  the  first  step  in  determining  the  anisotropy 
coefficients;  in  this  case,  however,  the  contravariant  velocities  are  used,  and  the 
derivatives  are  taken  with  respect  to  the  transformed  coordinates.  Also — and  this 
is  important — the  subgrid-scale  velocities  produced  by  the  turbulence  model  in 
this  case  are  the  contravariant  velocities. 

In  order  to  account  for  grid-generated  as  well  as  physical  anisotropies,  there¬ 
fore,  we  define  a  normalized  vector  of  anisotropy  coefficients,  C)  given  by  the 
formula: 


_  s\/3 

where  s  is  defined  according  to: 


(24) 


s 


t  = 


V3 


l|V(«c-i)|| 
llVUell  ■ 


(25) 


In  these  last  two  relations,  J~^  is  the  inverse  of  the  coordinate  transformation 
Jacobian  matrix,  i  is  a  unit  vector  and  iic  is  the  large-scale  contravariant  velocity 
vector,  filtered  to  preserve  only  the  smallest  (high-wavenumber)  components.  All 
norms  in  the  above  expressions  are  2-norms.  The  anisotropy  vector  is  normalized 
such  that  II  =  3;  thus  for  isotropic  turbulence  on  a  cubical  grid,  C '  *  =  1- 


5.3  Maps 

As  we  have  mentioned,  the  form  of  the  chaotic  maps  should  be  related  somehow 
to  the  statistical  properties  of  real  turbulence.  In  particular,  we  are  interested 
in  the  probability  density  of  the  velocity  fluctuations  Aw.  This  has  been  studied 
by  Kailasnath,  Sreenivasan,  k  Stolovitzky  [4],  in  research  that  was  also  funded 
by  the  AFOSR.  Led  by  Stolovitzky  [7]  (in  work  funded  by  DARPA),  these  same 
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authors  have  also  examined  Kolmogorov’s  refined  similarity  hypotheses.  The  first 
of  these  hypotheses  takes  the  form: 

Au{h)  =  (26) 

where  Ch,  is  the  local  average  of  the  energy  dissipation  rate,  over  an  interval  of  size 
h,  with  h  L.  Under  this  hypothesis,  U  is  a  stochastic  variable  with  a  probability 
density  function  (PDF)  that  depends  only  on  the  small-scale  Reynolds  number 
Reh  (Stolovitzky  et  al.use  a  slightly  different  Reynolds  number,  Rcr  =  Rei^^^). 
According  to  the  second  refined  hypothesis,  the  PDF  of  V  becomes  independent 
of  Rth,  and  thus  universal,  if  Reh.  »  1  (or  Rcr  >  1).  In  the  inertial  subrange,  the 
PDF  of  V  proves  to  be  nearly  Gaussian,  but  slightly  skewed,  unless  the  largest 
scales  are  excluded,  in  which  case  the  skewness  is  zero.  Juneja  et  o/.[3]  show  that 
the  variance  (V^)  should  be  approximately  2,  and  in  fact  Stolovitzky  et  al.[l] 
found  (K^)  =  1.88. 

Neglecting  anisotropy  effects,  our  subgrid-scale  model  for  the  velocity  generates 
values 

u,w(he)5M  (27) 

(where  i  indexes  the  grid  position)  so  there  should  be  a  close  relation  between  our 
stochastic  variable  M  and  the  stochastic  variable  V  appearing  in  Kolmogorov’s 
hypothesis.  In  fact,  the  velocity  increments  generated  by  our  model  take  the  form: 

(ui+i  -  u,_i)  =  C(/ie)i(Mi+i  -  M,_i),  (28) 

(with  C  =  0(1)).  The  mean  square  of  this  is 

((u.+i  -  u._i)')  =  C^hc)i{Ml,  -  2M.+iM._i  +  Ml,),  (29) 

but  if  Mt+i  and  Mi^i  are  independent  of  each  other,  the  mean  square  reduces  to 

((u.+i  -  u._i)')  =  C\he)^2Ml.  (30) 

Thus  the  variance  of  the  stochastic  variable  V  should  be  (approximately)  twice 
the  variance  of  our  stochastic  variable  M,  but  otherwise  the  two  variables  should 
share  the  same  PDF.  In  other  words,  the  PDF  of  the  stochastic  variable  M  should 
be  approximately  a  Gaussian,  with  unit  variance. 

Turbulent  velocity  fields  also  have  the  property  of  intermittent  small-scale  en¬ 
ergy  dissipation.  Juneja  et  a/.[3]  have  developed  a  technique  (so  far,  only  in  one 
dimension)  for  synthesizing  velocity  fields  that  have  many  of  the  same  statistical 
properties  as  physical  turbulence,  including  the  shape  of  the  PDF  for  the  veloc¬ 
ity  increments,  their  variance  and  skewness,  and  intermittent  dissipation.  Their 
procedure  is  based  on  a  multifractal  construction  involving  a  pair  of  unequal  mul¬ 
tipliers,  Cciscaded  through  several  stages.  Our  procedure  for  constructing  a  three- 
dimensional,  time- varying  velocity  field  on  the  small  scale  is  quite  different,  but 
we  have  suceeded  in  including  most  of  these  statistical  properties  in  our  model. 
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This  we  have  accomplished  by  constructing  our  stochastic  variable  M  from 
the  weighted  sum  of  chaotic  maps,  m^: 

M  =  Y,ajmj.  (31) 

3 

A  total  of  2^  maps  are  used  to  mimic  the  intermittent  dissipation  generated  by 
N  stages  in  the  construction  of  Juneja  tt  al..  Currently  we  are  using  four  maps 
(a  second-stage  construction);  in  this  case  the  appropriate  weights  aj  are: 


oi  =  1.2124,  02  =  0.7937,  03  =  0.5196,  04  =  0.7937.  (32) 

Just  one  of  these  maps  is  evaluated  at  each  iteration  of  the  stochastic  variable 
M.  At  each  grid-point  we  randomly  select  one  of  the  four  maps  to  begin  with, 
and  then  start  iterating  in  cyclic  sequence.  The  choice  of  which  map  to  evaluate 
at  any  given  iteration  could  also  (and  perhaps  should)  be  made  randomly,  but  at 
the  expense  of  slightly  greater  arithmetic  and  complexity,  which  at  present  we  are 
not  sure  is  justified.  In  any  ca^e,  the  procedure  appears  to  work  well.  For  a  case 
corresponding  to  fully  developed  turbulence — with  the  map  bifurcation  parameter, 
R  (see  below)  at  its  maximum  magnitude — Figure  2  shows  a  time-series  generated 
by  our  stochastic  variable  M,  while  Figure  3  shows  the  corresponding  probability 
density  function.  The  latter  quite  clearly  has  the  requisite  shape  and  variance. 

The  elemental  chaotic  maps  ruj  are  instances  of  the  normalized  quadratic  map: 


mi""'  =  H2{1  +  [l  -  \mf\  (l  +  i(l  +  ^/2)) 


(33) 
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or  the  normalized  tent  map: 


R{-2  -  3m5")) 
R{2  -  Zmf^) 


if  <  -j 
if  -  i  <  mf)  <  i 


(34) 


Here  (n)  is  the  iteration  number  and  R  is  the  map’s  bifurcation  parameter  (nor¬ 
malized  to  lie  between  ±1). 

In  order  to  relate  the  dynamics  of  the  chaotic  maps  to  the  dynamics  of  the  flow, 
we  make  the  bifurcation  parameter  R  a  function  of  a  local  large-scale  Reynolds 
number 


Re,  . 

1/ 


(35) 


where  I  is  a  characteristic  length  for  the  large-scale  motion.  This  Reynolds  num¬ 
ber  is  related  to  an  estimate  we  have  derived  for  the  ratio  of  the  Taylor  and 
Kolmogorov  microscales: 


A 

»7 


(36) 


Although  there  is  not  a  unique  expression  relating  R  and  Rcl,  the  following  ex¬ 
pression  behaves  reasonably  at  the  limits  and  matches  certain  critical  values  of 
the  Reynolds  number  and  bifurcation  parameter: 


R  “•  Rmax  tanh 


Rc  'll 

[\RecJ 

\  Umax  ^  _ 

(37) 
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In  this  expression,  Rtc  is  the  critical  value  of  Rei,  Rc  and  Rmax  are,  respec¬ 
tively,  the  critical  and  maximum  values  of  the  bifurcation  parameter  /?,  and  7 
is  a  constant  used  to  obtain  an  appropriate  transformation  from  physical  bifur¬ 
cation  parameters  to  those  of  the  chaotic  map.  In  one  of  our  present  studies, 
these  various  parameters  have  been  assigned  the  following  values:  Rec  «  1430, 
Rc  —  —0.2071,  Rmax  =  1,  and  7  =  3.6. 


5.4  Velocity  Correlations 

So  far  we  have  said  little  about  the  relative  independence  of  the  various  instances 
of  the  stochastic  variable,  M,  used  in  this  model.  At  present,  we  are  treating 
the  density  and  energy  fluctuations  as  if  they  were  independent,  although  they 
could  perhaps  be  related  through  the  energy  equation,  or  through  the  equation 
of  state.  In  fact,  we  anticipate  that  it  may  be  desirable  to  use  a  different  type  of 
chaotic  map  for  each  of  these  quantities,  each  map  being  tailored  to  best  match  the 
physical  behaviour  of  the  particular  quantity  being  modeled.  However,  we  know 
from  experiment  that  the  three  components  of  the  velocity  often  do  not  fluctuate 
independently.  This  interdependence  can  be  modeled  by  coupling  together  the 
three  maps  used  to  generate  the  velocity  fluctuations. 

For  now,  let  us  consider  the  correlations  between  the  small-scale  velocity  com¬ 
ponents  u'f.  at  a  single  point  in  space.  These  form  a  tensor: 


Rkt  = 


(38) 


where  each  of  the  elements  has  been  normalized  to  lie  in  [— 1,-fl].  Once  again 
considering  velocity  fluctuations  on  a  scale  just  above  the  grid  scale,  we  estimate 


(Vufc)  •  (Vuf) 

||Vu,||  ||Vu,||- 


(39) 


(As  in  the  section  5.2,  these  gradients  will  be  evaluated  in  terms  of  the  filtered 
contravariant  velocities  and  the  transformed  coordinates.)  The  problem  is  to  de¬ 
fine  a  coupling  between  the  stochastic  variables  for  the  three  velocity  components 
that  will  reproduce  this  tensor. 

Recall  from  equation  (10)  that  the  subgrid-scale  velocity  vector  Ws  (or  us  in 
the  incompressible  case)  depends  on  a  vector  of  stochastic  variables  M.  Let  this 
latter  vector  have  components 


M  =  [Afi,Af2,A^3f.  (40) 

Furthermore,  let  each  of  the  components  Mk  be  composed  of  the  weighted  sum 
of  three  instances  of  the  stochastic  variable  M: 


A4/;  =  ockiMi  -f  Oik2^2  +  QJJbsAfa, 


(41) 
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where  Mi,  M2,  and  M3  are  three  statistically  independent  instances  of  the  stochas¬ 
tic  variable  defined  by  equation  (31);  that  is 


{MmMn) 


(M^)  if  m  =  n, 
0  if  m  n. 


And  let  us  normalize  the  coefficients  so  that 


(42) 


E“L  =  1-  (43) 

m 

Then  if  the  correlation  tensor  Rkt  is  formed  from  the  modelled  small-scale  veloc¬ 
ities,  the  result  is 

=  =  (44) 


Equating  this  with  (39)  is  accomplished  straightforwardly  by  setting 


^km  — 


llViill  • 


(45) 


This  works  because  the  interaction  between  the  averaging  process,  the  statistical 
independence  of  the  maps  Mm,  and  the  form  of  (41)  mimics  the  dot  product 
appearing  in  (39). 

Note  that,  in  a  turbulent  flow,  the  components  of  Vu  are  likely  to  fluctuate. 
Consequently,  the  coupling  coefficients  a,*  will  fluctuate  also.  Under  these  cir¬ 
cumstances,  the  coupling  coefficients  will  still  contribute  to  the  correlation  tensor 
to  the  extent  that  they  really  are  correlated. 


5.5  Time-Scale 

In  order  to  use  our  subgrid-scale  model  in  a  computational  simulation  of  turbu¬ 
lence,  we  also  must  relate  the  iterations  of  the  stochastic  variable  to  the  time-step 
k  used  in  the  simulation,  and  to  a  characteristic  time  for  the  subgrid-scale  eddies. 
The  latter  we  estimate  ais  a  turnover  time  for  a  circular  eddy  of  diameter  hk'. 

where  Uk  is  the  circumferential  velocity  of  the  eddy  and  /m  is  a  fundamental 
freqency  characteristic  of  the  chaotic  maps  that  make  up  the  stochastic  variable 
M. 

Now  in  the  isotropic  case, 

hk  =  (3)-U,  (47) 
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while 


\ 


(48) 


Therefore,  after  some  manipulation,  we  derive 

_  tt/mRcI 

y/SC4Vu\\' 


(49) 


as  an  estimate  for  the  subgrid  time-scale.  This  time-scale  can  also  be  estimated  in 
terms  of  the  grid-scale  vorticity,  with  results  that  are  comparable  in  the  inertial 
subrange  of  turbulence. 

The  subgrid  time-scale  tg  proves  to  be  somewhat  larger  than  the  time-steps 
that  typically  must  be  used  in  the  numerical  simulation  of  fluid  flows.  Thus  one 
iteration  of  the  stochastic  variable  M  corresponds  to 


tl  - 

k  ~  kV3Cu\\Vu\\ 


(50) 


time  steps.  The  singularity  that  would  occur  when  l|Vw||  =  0  is  not  a  problem, 
because  the  model  is  turned  off  when  Rth  <  Rckmin  —  <^(1)-  Practically,  under 
this  circumstance,  one  limits  n  to  being  less  than  some  maximum  value,  and 
rechecks  llVw||  at  each  time-step. 


5.6  Mass  Conservation 

For  incompressible  flows,  a  final  step  in  the  construction  of  the  subgrid-scale  ve¬ 
locity  field  is  the  removal  of  the  divergence,  i.e.  the  projection  of  this  velocity 
onto  a  divergence-free  subspace  of  three-dimensional  vector  fields.  This  is  justi¬ 
fied  heuristically  on  the  grounds  that  the  difference  between  the  longitudinal  and 
transverse  velocity  correlations  in  a  turbulent  flow  is  affected  by  the  divergence- 
free  constraint  (see,  for  example,  Batchelor  [1]).  We  have  also  found  that  our 
computational  results  seem  better  when  we  take  this  step.  On  this  basis  we  also 
believe  that  mass  conservation  should  be  imposed  on  the  small-scale  field  in  com¬ 
pressible  flow.  This  remains  an  area  we  are  actively  researching. 

5.7  Summary  of  the  Subgrid-Scale  Model 

In  summary,  we  mimic  subgrid-scale  turbulence  by  means  of  stochastic  small-scale 
fields,  each  of  which  is  produced  by  generating,  at  each  point  on  a  grid,  a  value 
that  is  the  product  of  amplitude  factors  with  normalized  stochastic  variables.  The 
basic  amplitude  factors  are  determined  in  accordance  with  Kolmogorov’s  theory 
of  turbulence,  although  we  have  been  forced  to  rely  on  some  heuristic  arguments 
in  order  to  estimate  the  amplitude  of  the  density  and  energy  fluctuations.  Our 
stochastic  variables  are  constructed  from  a  combination  of  chaotic  maps,  and 
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their  statistics  are  tailored  to  match  at  least  the  basic  statistics  of  physical  tur¬ 
bulence  in  the  inertial  subrange.  Anisotropy  in  the  velocity  field,  as  well  as  the 
cross-correlation  of  the  velocity  field  components,  are  introduced  as  appropriate 
according  to  a  local  sample  of  the  structure  of  the  large-scale  velocity  field.  A 
time-scale  is  determined  (locally)  for  the  subgrid-scale  flow;  this  time-scale  deter¬ 
mines  in  turn  how  often  the  stochastic  fields  are  updated.  Finally,  care  is  taken 
to  ensure  that  the  subgrid-scale  flow  is  mass-conserving. 

6  Implementation  for  Compressible  Flow 

The  first  implementation  of  the  new  turbulence  modeling  scheme  for  compressible 
flow  began  in  May  1993.  Our  goal  was  to  construct  the  implementation  as  a 
relatively  simple  modification  to  the  existing  finite  volume  compressible  flow  solver 
GASP  and  to  test  the  scheme  on  a  supersonic  compression  ramp  flow.  This  was 
carried  out,  based  on  the  scheme  as  described  in  the  Mechanical  Engineering 
Department  CFD  Report  by  Hylin  and  McDonough  in  April  1994.  Although  the 
small-scale  fluctuations  appeared  to  behave  exactly  as  designed,  no  large-scale 
turbulent  boundary  layer  structures  appeared  in  the  computed  solution.  It  is 
crucial  to  represent  such  structures  in  the  large-scale  solution  because  they  are 
responsible  for  the  majority  of  the  turbulent  transport  normal  to  the  wall  in 
boundary  layer  flows.  It  is  not  clear  whether  our  small-scale  modeled  quantities 
must  actually  generate  the  large-scale  structures  or  merely  react  to  them,  but  there 
is  strong  evidence  in  the  literature  that  small-scale  coherent  motions  within  100  y'*' 
units  of  the  wall  are  responsible  for  the  maintenance  of  the  large-scale  outer-region 
turbulence.  We  performed  numerical  experiments  in  which  larger-than-grid-size 
versions  of  these  near-wall  coherent  motions  were  imposed  as  perturbations  to 
initially  laminar  supersonic  flat  plate  and  compression  ramp  flows.  In  both  cases, 
the  computed  solutions  contained  large-scale  structures  of  the  types  observed  in 
physical  turbulent  flows;  this  is  an  example  of  effective  interaction  of  a  large- 
scale  solution  and  a  smaller-scale  perturbation,  a  model  for  our  scheme.  We 
concluded  that  it  is  indeed  possible  to  induce  large-scale  turbulent  structures  in 
the  large-scale  solution,  and  that  for  this  to  occur,  the  small-scale  fields  must 
contain  spatially  coherent  structures. 

Based  on  these  conclusions  and  continuing  developments  in  the  derivation  of 
the  underlying  scheme,  an  improved  implementation  is  in  preparation  at  present. 
The  new  implementation  is  much  simpler  to  retrofit  to  existing  flow  solvers  and 
requires  much  less  computer  memory  to  run  than  the  first  implementation.  Con¬ 
servation  of  mass  will  be  imposed  on  the  small-scale  fields,  and  at  least  some 
calculations  will  include  upstream  boundary  conditions  designed  to  induce  the 
large-scale  turbulent  structure  in  the  solution. 

In  following  sections  we  will  describe  the  work  to  date. 
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Figure  4:  Instantaneous  complete  density  field. 

7  First  Implementation 

In  the  first  months  of  this  effort,  the  details  of  the  original  modeling  scheme  were 
developed  and  refined  at  the  same  time  we  studied  the  GASP  source  code  in  prepa¬ 
ration  for  incorporation  of  the  new  scheme.  Based  on  our  initial  understanding 
of  the  function  of  the  small-scale  fluctuations  in  the  method,  distinct  small-scale 
variables  were  assigned  to  cell  centers  and  cell  faces  to  provide  for  variation  of  the 
small-scale  solution  across  individual  cells.  We  formulated  an  acceptable  method 
of  adapting  the  Roe  flux  difference  scheme  available  in  GASP  to  solution  of  the 
ATD  large-scale  equations  as  they  were  formulated  at  the  time.  This  was  incor¬ 
porated  into  the  code.  New  subroutines  to  generate  the  small-scale  quantities  a.s 
well  as  the  auxiliary  variables  in  terms  of  generalized  coordinates  were  written. 
We  created  new  code  input  variables  and  source  code  to  control  the  operation 
of  the  small-scale  model  and  the  production  of  intermediate  visualization  output 
and  restart  files  for  large-scale  and  small-scale  variables. 

The  new  scheme  was  tested  on  a  2-D  compression  ramp  flow.  The  correspond¬ 
ing  physical  experiments  were  conducted  at  Princeton.  The  geometry  is  a  24° 
ramp  in  a  Mach  2.85  flow  at  a  Reynolds  number  of  1.6  x  10® /in.  The  calculations 
were  carried  out  on  a  240x80  grid  that  was  highly  stretched  in  the  wall-normal 
direction.  For  purposes  of  reference  in  the  following  figures,  the  ramp  corner  is 
located  at  cell  89  in  the  streamwise  direction.  Figure  4  shows  an  instantaneous 
complete  density  field  during  the  calculation;  this  shows  the  main  shock  structures 
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Figure  5;  Instantaneous  velocity  fields  in  compression  corner;  (a),  complete  ve¬ 
locity;  (b),  large-scale  velocity;  (c),  small-scale  velocity. 


and  the  recirculation  region.  The  shock  structure  fluctuates  in  the  streamwise  di¬ 
rection  in  the  experiments,  but  no  such  motion  was  evident  in  our  calculations. 

The  small-scale  fluctuations  behaved  as  designed,  with  appropriate  dynamics 
and  amplitudes  determined  by  local  large-scale  conditions. 

Figure  5  shows  a  close-up  view  of  the  compression  corner,  with  every  fifth 
velocity  vector  in  the  complete,  large-scale,  and  small-scale  velocity  fields;  the 
recirculation  region  is  evident  in  (a),  the  complete  velocity  field.  Although  it  is 
not  obvious,  somewhat  less  grid-scale  oscillations  are  present  in  the  large-scale 
field  (b)  than  in  the  complete  field.  Panel  (c)  shows  the  corresponding  small-scale 
velocity  field  and  its  contribution  to  the  complete  field.  The  vectors  in  this  panel 
have  been  scaled  up  by  a  factor  of  about  two  from  that  used  in  the  other  panels  to 
make  them  more  visible.  For  the  small-scale  field,  notice  that  the  largest  amplitude 
fluctuations  are  located  close  to  the  wall  where  the  shear  is  highest.  This  feature 
is  built  into  the  model  for  the  velocity  fluctuations.  Figure  6  illustrates  this  again; 
in  the  figure,  the  vectors  have  been  scaled  up  by  a  factor  of  five  from  the  previous 
figure  to  make  them  more  visible.  Note  that  the  modeled  small-scale  velocity 
fluctuations  are  essentially  nonexistent  in  the  freestream  regions  on  the  flow  field 
as  they  should  be. 

The  next  figure  shows  line  plots  of  instantaneous  complete  primitive  variable 
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Figure  6:  Instantaneous  small-scale  velocity  field  (scaled  lOx),  full  domain. 

quantities  along  a  streamwise  line  located  at  the  fifth  cell  from  the  wall.  The 
fluctuation  amplitudes  vary  widely  along  the  domain,  but  they  are  of  reasonable 
size  in  comparison  with  local  complete  values.  The  evident  spatial  chaos  is  ac¬ 
companied  by  temporal  chaos,  although  no  time  series  are  presented  here. 

Results  from  these  calculations  were  presented  at  meetings  of  the  American 
Physical  Society  Division  of  Fluid  Dynamics  in  Albuquerque,  NM,  in  November 
1993,  and  in  Atlanta,  GA,  in  November  1994. 

8  Induction  of  Large-Scale  Turbulent  Boundary 
Layer  Structure 

As  promising  as  the  model  behavior  was  in  the  ramp  calculations,  the  interaction 
of  the  fluctuations  with  the  large-scale  solution  was  less  than  satisfactory.  No 
large-scale  coherent  boundary  layer  structures  appeared  in  the  time-dependent 
ATD  solution.  This  was  considered  a  serious  defect  since  such  structures  are 
responsible  for  the  majority  of  the  turbulent  transport  and  they  are  the  main 
features  to  be  computed  in  the  large-scale  solution.  No  statistics  were  computed 
for  the  2-D  ramp  flowfield  because  the  absence  of  large-scale  motion  in  the  field 
removes  any  hope  of  a  reasonable  comparison  with  experiment. 

Attention  focussed  on  how  to  ensure  proper  interaction  of  the  small-scale  quan¬ 
tities  with  the  large-scale  solution  in  the  context  of  a  finite  volume  scheme.  The 
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Figure  8:  (a),  grid;  (b),  boundary  conditions. 


consensus  of  the  coherent  structure  literature  and  recent  developments  in  stabil¬ 
ity  theory  is  that  interaction  of  quasi-streamwise  vortices  and  streaks  in  the  inner 
near-wall  layer  are  likely  responsible  for  initiating  and  maintaining  wall-bounded 
turbulence,  and  in  particular,  this  interaction  directly  produces  the  large-scale 
structures  that  we  must  compute.  We  carried  out  3-D  computations  of  a  super¬ 
sonic  flat  plate  boundary  layer  under  the  same  flow  conditions  used  in  the  ramp 
calculation. 

The  computational  grid  contained  120  x  60  x  20  cells.  As  shown  in  Figure 
8,  the  laminar  flow  initial  condition  was  perturbed  by  streamwise  vortices  im¬ 
posed  as  a  stand-in  for  the  small-scale  quantities  we  seek  to  model.  No  turbulence 
model  was  employed  in  the  calculation.  The  experiments  showed  that  such  a  com¬ 
puted  boundary  layer  is  indeed  susceptible  to  “transition”  under  the  influence  of 
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the  streak-vortex  interaction  mechanism.  The  large-scale  coherent  structures  in 
the  solution  were  quite  similar  in  character  to  those  observed  in  fluid  turbulence. 
Figure  9(a),  a  contour  plot  of  streamwise  momentum  at  the  first  cell  center  off 
the  plate  surface  (plane  j=l),  thresholded  to  highlight  values  less  than  20%  of  the 
freestream  value,  shows  the  streaky  structures  described  as  essentially  ubiquitous 
in  the  sublayer  and  buffer  regions  of  all  turbulent  boundary  layers.  These  struc¬ 
tures  occur  at  very  small  length  scales  near  the  wall  in  physical  flows,  but  they 
appear  at  scales  larger  than  the  grid  size  in  the  calculation.  Figure  9(b)  is  a  plot 
of  wall-pressure  contours,  limited  to  those  greater  than  the  initial  freestream  pres¬ 
sure.  The  form  of  these  contours  corresponds  to  literature  descriptions  of  the  high 
wall-pressure  regions  as  relatively  rounded  spots  rather  than  elongated  regions. 

Figure  10  shows  contours  of  spanwise  vorticity  in  the  wall-normal  plane  on  the 
domain  spanwise  centerline.  The  four  parts  of  the  figure  are  from  successive  times 
in  the  flow  evolution.  The  most  noticeable  feature  of  this  sequence  is  the  devel¬ 
opment  and  convection  downstream  of  a  large  bulge  in  the  rotational-irrotational 
interface,  and  an  accompanying  incursion  of  irrotational  fluid  toward  the  wall. 
This  is  a  common  feature  often  observed  in  turbulent  boundary  layers.  The  con¬ 
vection  speed  of  the  bulge  can  be  estimated  from  this  sequence  to  be  about  80% 
of  the  freestream  velocity;  this  is  consistent  with  values  found  in  incompressible 
low-Reynolds  number  boundary  layers,  but  is  lower  than  the  corresponding  per¬ 
centage  for  boundary  layers  at  this  Mach  number,  which  has  been  measured  to  be 
about  90%.  In  the  figure  we  also  observe  sloping  regions  of  very  high  shear  near 
the  wall.  These  shear  layers  also  agree  well  with  experimental  observations. 

Figure  11(a)  is  similar  to  Figure  10(c);  it  shows  spanwise  vorticity,  with  fea¬ 
tures  evident  as  before.  We  call  attention  to  some  details  here.  Figure  11(b)  shows 
velocity  vectors  in  the  plane.  The  profiles  at  any  given  downstream  station  bear 
little  resemblance  to  the  nearly  undisturbed  profiles  at  the  upstream  boundary  at 
left.  In  the  downstream  third  of  the  domain,  we  find  a  structure  with  a  sharply- 
defined  sloping  upstream  interface  that  appears  to  be  almost  a  discontinuity  in 
the  streamwise  velocity.  The  structure  extends  all  the  way  across  the  boundary 
layer.  Regions  of  both  high-  arid  low-speed  fluid  exist  within  this  structure,  with 
a  range  of  scales  present.  This  corresponds  to  the  feature  noted  by  many  to  be  the 
most  readily  detected  in  experiments,  the  sloping  J-scale  “backs”  of  large  bulges. 
The  fluid  ahead  of  such  backs  is  frequently  observed  to  rotate  slowly  in  a  direction 
consistent  with  its  downstream  motion.  In  Figure  11(c)  we  plot  the  same  velocity 
vector  field,  but  here  in  a  reference  frame  moving  at  the  previously  estimated 
convection  velocity  of  the  structure.  The  field  does  in  fact  show  slow  rotation  as  if 
rolling  downstream  along  the  plate.  Further,  at  two  locations  along  the  upstream 
side  of  these  structures,  high-speed  fluid  from  upstream  encounters  low-speed  fluid 
with  the  consequence  of  very  high  local  pressure.  This  is  illustrated  for  this  plane 
in  Figure  11(d),  which  shows  contours  of  pressure.  The  two  dark  regions  in  the 
plot  about  halfway  up  in  the  domain  are  the  areas  in  question. 

Figure  12(a)  shows  cross-stream  sections  at  the  five  streamwise  stations  marked 
in  Figure  11.  Figure  12(b)  is  a  cross-stream  section  of  the  entire  domain  at  a:  =  3.5 
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Figure  9:  Wall  contours  at  t  =  0.50  msec:  (a),  low  streamwise  momentum;  (b), 
high  pressure. 


(d) 

Figure  10:  Contours  of  spanwise  vorticity  at  time  =  (a),  0.40  msec;  (b),  0.45  msec; 
(c),  0.50  msec;  (d),  0.55  msec. 

in,  showing  contours  of  pressure  and  projections  of  the  velocity  vectors  onto  this 
plane.  The  initial  streamwise  vortices  are  still  very  much  in  evidence,  but  smaller 
scale  motions  have  also  developed.  Note  that  the  pressure  contours  in  this  plane 
show  the  existence  of  low-pressure  cores  in  the  vortices;  the  low-pressure  cores 
have  been  observed  to  be  a  reliable  and  convenient  indicator  of  vortical  structures 
in  complex  flowfields.  Vortical  structures  are  connected  in  some  way  to  each  of  the 
other  classes  of  motion;  in  the  viscous  sublayer  and  buffer  layer,  quasi-streamwise 
vortices  are  observed  to  play  a  central  role,  although  the  precise  dynamics  have 
not  yet  been  established.  In  the  sections  of  Figure  12(a)  we  see  just  such  a 
vortical  structure.  The  sections  show  velocity  vectors  projected  onto  the  cross¬ 
stream  plane,  with  pressure  contours  superimposed.  In  each  section  we  see  flow 
circulating  around  a  core  of  low  pressure.  The  vortex  core  is  located  within  two 
or  three  grid  cells  of  the  wall  in  the  upstream  section  at  left.  The  core  rises  going 
downstream,  in  accordance  with  the  observation  that  quasi-streamwise  vortices  are 
quite  small  and  nearly  horizontal  near  the  wall,  but  their  diameter  and  inclination 
increase  downstream,  until  in  the  mid-ranges  of  the  layer,  the  inclination  angle 
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(d) 

Figure  11:  Time  =  0.50  msec;  (a),  spanwise  vorticity  contours;  (b),  velocity  vec¬ 
tors;  (c),  velocity  vectors  in  reference  frame  convecting  at  80%  Uoo‘,  (d),  pressure 
contours. 
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Figure  12:  Time  0.50  msec;  (a),  sections  through  streamwise  vortex  at  five  stream- 
wise  stations;  cross-stream  velocity  vectors  and  pressure  contours;  (b),  cross¬ 
stream  velocity  vectors  and  pressure  contours  at  streamwise  station  x  =  3.5  in. 
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averages  about  45°-60°. 

The  flowfield  appears  to  have  some  structural  features  that  occur  within  the 
viscous  sublayer  and  buffer  layer.  These  structures  appear  at  scales  much  larger 
than  those  expected  in  the  physical  flow.  The  presence  of  these  structures  in 
our  computed  solution  is  consistent  with  conjectures  about  the  importance  of 
interaction  of  vortices  and  streaks  in  physical  wall-bounded  turbulence.  Such  an 
interaction  will  likely  be  necessary  in  the  large-scale  computed  solution  for  large- 
scale  ‘turbulent’  motion  to  be  maintained,  unless  the  modeled  turbulent  sub-grid 
fluctuations  can  be  made  to  carry  out  this  function  in  some  way. 

On  the  other  hand,  the  coherent  structures  that  typically  occur  at  (J-scale,  such 
as  the  large  bulges  and  the  sloping  “backs”,  were  well-represented  in  size  and  speed 
of  propagation.  This  suggests  that  we  may  hope  to  compute  such  structures,  our 
ultimate  goal,  in  the  future  with  an  appropriate  numerical  scheme. 

The  turbulent  structures  in  this  computed  solution  were  provoked  by  the  pres¬ 
ence  of  perturbations  in  the  form  of  streamwise  vortices.  We  infer  that  if  the  mod¬ 
eled  small-scale  variable  fields,  which  function  similarly  to  these  perturbations  in 
our  turbulence  modeling  scheme,  contain  some  component  of  such  motion,  that  a 
similar  transition  effect  may  occur.  This  is  one  of  the  bases  for  our  belief  that  co¬ 
herent  structure  will  be  necessary  in  the  small-scale  fields  for  effective  interaction 
with  the  large-scale  solution. 

The  results  of  these  experiments  were  reported  in  more  detail  in  an  AIAA 
Paper  presented  at  the  34th  AIAA  Aerospace  Sciences  Meeting,  Reno,  NV,  in 
January  1996. 

A  similar  experiment  was  carried  out  for  the  3-D  compression  ramp  flow  that 
was  used  as  our  first  test  case.  In  this  calculation,  the  coiriputational  grid  con¬ 
tained  160  X  50  X  20  cells.  As  in  the  previous  experiment,  streamwise  vortex 
perturbations  were  superimposed  on  the  laminar  supersonic  initial  condition  and 
the  flow  was  evolved  without  a  turbulence  model.  Turbulent  structure  developed 
quickly  and  did  not  seem  to  diminish  as  the  calculation  proceeded,  in  contrast 
with  the  flat  plate  case.  Activity  downstream  of  the  recirculation  shock  was 
particularly  dramatic.  Animation  of  the  solution  revealed  time-dependent  shock 
structures  that  compared  well  with  the  visual  observations  of  experimentalists. 
A  typical  field  is  presented  in  Figure  13j  the  figure  shows  filled  density  contours 
at  an  instant  about  half  a  millisecond  (physical  time)  into  the  calculation.  The 
large-scale  vortical  structures  rolling  up  in  the  boundary  layer  downstream  of  the 
recirculation  shock  are  filled  with  low-speed  fluid  that  produces  secondary  shocks 
as  high-speed  fluid  strikes  them  from  behind.  One  highly  unsatisfactory  aspect  of 
the  computed  solution  was  that  the  recirculation  shock  moved  upstream  during 
the  calculation  until  it  almost  reached  the  inflow  boundary.  The  more  or  less  final 
location  of  this  shock  was  very  far  from  that  measured  in  the  Princeton  lab. 

The  results  of  this  calculation  were  reported  at  the  48th  Annual  Meeting  of  the 
American  Physical  Society  Division  of  Fluid  Dynamics,  Irvine,  CA,  in  November 
1995. 

The  results  of  this  calculation  were  presented  at 
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Figure  13:  Instantaneous  density  field,  24°  compression  ramp. 

9  Revised  Implementation 

At  present,  revisions  to  the  formulation  and  implementation  of  the  turbulence 
modeling  scheme  are  in  progress.  The  underlying  formulation  of  the  decomposed 
equations  has  been  reconsidered.  The  formulation  for  the  large-scale  equations 
now  requires  construction  of  the  large-scale  part  of  the  advective  terms  in  the 
governing  equations,  with  an  approximation  to  the  complete  solution  used  in  this 
process.  This  removes  the  question  of  how  to  split  the  various  scale  interaction 
terms  among  the  several  sets  of  decomposed  equations  and  avoids  the  accompany¬ 
ing  theoretical  difficulties  of  deriving  characteristic-based  solution  schemes  with 
rigor.  From  the  practical  point  of  view,  it  also  greatly  simplifies  implementation; 
the  numerical  methods  in  the  existing  flow  solver  can  be  used  for  the  new  formu¬ 
lation  without  modification  except  that  the  large-scale  and  small-scale  quantities 
must  be  combined  to  produce  approximate  complete  quantities  at  the  beginning 
of  each  time  step. 

Our  computed  compressible  flow  experiments  with  streamwise  vortex  pertur¬ 
bations  and  our  experience  with  calculating  incompressible  flow  over  a  backward¬ 
facing  step  lead  us  to  believe  that  effective  interaction  between  the  modeled  small- 
scale  quantities  and  the  computed  large-scale  solution  occurs  when  the  small-scale 
field  contains  structures  at  scales  larger  than  the  grid  scale.  This  implies  that  it 
is  not  necessary  to  maintain  small-scale  variables  at  cell  centers  and  cell  faces  as 
we  have  done  previously,  because  variation  across  individual  cells  has  little  effect 
on  the  large-scale  solution.  In  addition,  use  of  small-scale  variables  at  several  lo¬ 
cations  per  cell  greatly  increases  the  memory  requirements  over  the  original  flow 
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solver.  Therefore,  in  the  new  implementation,  small-scale  variables  are  carried 
only  at  cell  centers,  the  nominal  locations  for  the  large-scale  variables. 

Our  concern  with  effective  scale  interaction  will  also  affect  our  imposition  of 
the  continuity  equation  on  the  modeled  small-scale  field.  We  will  prefer  a  method 
that  automatically  produces  coherent  structure  in  the  small-scale  velocity  field 
itself  at  each  time  step,  in  the  same  manner  as  in  the  incompressible  case. 

Experiments  with  the  reformulated  scheme  will  be  reported  in  future  publica¬ 
tions. 
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Weatherly,  D.  C.  and  McDonough,  J.  M.,  “Computed  Large-Scale  Compressible 
Boundary  Layer  Structure  Stimulated  by  Vortical  Perturbations,”  presented  at 
AIAA  34th  Aerospace  Sciences  Meeting,  Reno,  Jan.  15-18,  1996. 

Conference  Presentations 

Hylin,  E.  C.,  Weatherly,  D.  C.,  and  McDonough,  J.  M.,  “Modeling  the  Subgrid- 
Scale  Flow  with  a  Chaotic  Map,”  presented  at  Amer.  Phys.  Soc.  46th  Ann.  Mtg. 
Division  of  Fluid  Dynamics,  Albuquerque,  NM,  Nov.  21-23,  1993. 

Weatherly,  D.  C.,  Hylin,  E.  C.  and  McDonough,  J.  M.,  “Additive  Turbulent  De¬ 
composition  with  Chax)tic  Map  Subgrid-Scale  Model  for  Compressible  Flow,”  pre¬ 
sented  at  Amer.  Phys.  Soc.  46th  Ann.  Mtg.  Division  of  Fluid  Dynamics, 
Albuquerque,  NM,  Nov.  21-23,  1993. 

McDonough,  J.  M.,  Zhong,  X.  and  Hylin,  E.  C.,  “Additive  Decomposition  of 
the  Navier-Stokes  Equations  with  Chaotic  Map  Subgrid  Models:  Application  to 
Intermittent  Pipe  Flow,”  presented  at  Amer.  Phys.  Soc.  47th  Ann.  Mtg.  Division 
of  Fluid  Dynamics,  Atlanta,  GA,  Nov.  20-22,  1994. 

Weatherly,  D.  C.,  Hylin,  E.  C.  and  McDonough,  J.  M.,  “Additive  Turbulent  De¬ 
composition  of  the  Navier-Stokes  Equations  with  Chaotic  Map  Subgrid  Models: 
Application  to  Supersonic  Ramp  Flow,”  presented  at  Amer.  Phys.  Soc.  47th 
Ann.  Mtg.  Division  of  Fluid  Dynamics,  Atlanta,  GA,  Nov.  20-22,  1994. 

McDonough,  J.  M.,  Yang,  Y.  and  Hylin,  E.  C.,  “Incompressible  Flow  over  a 
Backward-Facing  Step,”  presented  at  Amer.  Phys.  Soc.  48th  Ann.  Mtg.  Di¬ 
vision  of  Fluid  Dynamics,  Irvine,  CA,  Nov.  19-21,  1995. 

Weatherly,  D.  C.,  Hylin,  E.  C.  and  McDonough,  J.  M.,  “Supersonic  Compres¬ 
sion  Ramp  Calculations  with  Dependent- Variable  Sub-grid  Models,”  presented  at 
Amer.  Phys.  Soc.  48th  Ann.  Mtg.  Division  of  Fluid  Dynamics,  Irvine,  CA,  Nov. 
19-21,  1995. 
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Invited  Talks  and  Lectures 

McDonough,  J.  M.,  “The  Additive  Turbulent  Decomposition  Formalism:  An  Al¬ 
ternative  Approach  to  Turbulence  Modeling  and  Simulation,”  presented  at  Hong 
Kong  Univ.  Sci.  Tech.,  Clear  Water  Bay,  Hong  Kong,  May  24,  1994. 

_ ,  “Additive  Turbulent  Decomposition  of  the  Navier-Stokes  Equations:  An 

Alternative  Formalism  for  Turbulence  Modeling,”  presented  at  UCLA  Mechanical, 
Aerospace  k  Nuclear  Engineering  Seminar,  Mar.  23,  1995. 

_ ,  “Turbulence  Computations  Based  on  Unaveraged  Equations  and  Non¬ 
linear  Chaotic  Map  Subgrid-Scale  Models,”  presented  at  NASA  Lewis  Research 
Center,  Heat  Transfer  Branch,  Aug.  22,  1995. 

_ ,  “Small-Scale  Turbulence  Modeling  and  Simulation  via  Additive  Turbu¬ 
lence  Decompositions  of  the  Navier-Stokes  Equations,”  presented  at  the  Institute 
for  Mathematics  and  its  Applications,  Univ.  of  Minnesota,  Oct.  26,  1995. 
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